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1  Introduction 

The  knowledge  of  a  spacecraft’s  position  in  it’s  orbit  has  always  been  of 
paramount  importance,  and  the  estimation  of  the  orbit  has  historically  been 
a  ground-based  operation.  The  computation  facilities  required  to  handle  the 
sophisticated  optimal  estimation  algorithms  have  been  of  the  VAX  class,  or 
larger,  and.  of  necessity,  been  confined  to  well-ventilated  computer  centers. 
This  is  changing.  CPUs  which  can  outperform  a  VAX  11/780  by  factors 
of  ten  or  twenty  are  now  available  in  a  two  chip  set  which  consume  only 
a  few  watts.  Orbit  estimation,  now  being  called  Autonomous  Navigation 
(AutoNav).  can  now  be  physically  performed  on  its  associated  vehicle.  This 
paper  reports  on  our  investigations  of  one  such  AutoNav  scheme. 

AutoNav  systems  for  spacecraft  comprise  four  major  components: 

Reference  Points:  Spacecraft  orbit  estimation  requires  knowledge  of  refer¬ 
ence  vectors  to  targets  or  bodies  whose  positions  are  known  in  inertial 
space.  These  measurements  are  sometimes  used  in  conjunction  with 
the  knowledge  of  satellite  attitude  to  determine  the  orbital  parameters 
of  the  satellite.  These  known  landmarks  may  be  celestial  bodies  such 
as  the  sun.  moon,  or  planets.  They  may  be  known  points  on  the  Earth 
or  moon,  or  they  may  be  other  satellites  in  known  orbits,  such  as  GPS. 
Without  compromising  the  utility  of  the  general  technique  described, 
this  study  has  considered  the  use  of  Earth-based  landmarks  to  define 
a  system  that  is  easily  realisable. 

Sensor(s):  Depending  on  the  type  of  reference  points  used,  sensors  which 
observe  these  references  may  be  ranging  devices  (e.g.,  when  used  with 
GPS  or  cooperative  ground  based  transmitters)  or  angle  sensing  de¬ 
vices.  In  this  study,  the  latter  type  was  used;  an  optical  sensor  of  the 
NRL  Rad-Hard  Star  Tracker  class. 


Measurement  Processor:  The  measurements  made  by  the  sensor(s)  on 
the  reference  point(s)  must  be  converted  into  corrrections  to  the  Au- 


toNav  system’s  opinion  of  where  it  is.  The  ideal  algorithm  for  this  is 
the  Extended  Kalman  Filter. 

Orbit  Propagator:  When  the  reference  points  are  not  visible  to  the  Au- 
toNav  sensor,  and  measurements  are  not  available,  the  AutoNav  system 
must  be  able  to  project  its  position  as  time  passes  until  new  measure¬ 
ments  come  in.  This  is  the  job  of  the  Orbit  Propagator. 

These  four  components  are  described  much  more  fully  in  following  sections. 

There  are  two  approaches  to  landmark  (Earth  surface)  tracking: 

Imaging  Sensor:  This  concept  involves  the  use  of  an  on-board  imaging 
sensor,  such  as  SPOT  or  other  Earth  imaging  sensor,  to  extract  the 
position  of  known  topographic  features.  While  such  landmark  data  has 
been  used  to  improve  the  attitude  estimate  of  some  systems,  the  use 
of  this  method  has  not  been  demonstrated  in  an  operational  attitude 
control  system.  It  is  assumed  that  high  resolution  can  be  achieved  only 
after  ground  processing  of  the  imager  data.  This  limits  the  application 
of  this  type  of  data  in  a  real-time  autonomous  navigation  system. 

Beacon/Target  Sensor:  Tracking  with  ground-based  beacons  involves  us¬ 
ing  a  point  source  tracker  to  acquire  and  track  a  cooperative  target  or 
beacon  on  the  ground  at  a  known  location.  The  target  could  be  a  laser 
beacon  source  or  a  passive  retroreflector  used  with  an  on-board  illu¬ 
minator.  This  concept  can  be  extended  to  tracking  beacons  on  other 
sources,  such  as  the  moon,  aircraft,  ships,  or  other  satellites. 

Sensor  design,  operation,  and  accuracy  are  essentially  identical  to  those 
of  a  star  tracker,  allowing  application  of  a  well-developed  technology.  The 
beacon  could  employ  a  visible  or  near-infrared  laser  source,  such  as  a  laser 
diode,  for  use  with  a  tracker  employing  a  silicon  charge  transfer  device  (CTD ) 
detector.  The  use  of  an  active  sensor  illuminating  a  passive  retroreflector  on 
the  ground  has  advantages  in  autonomy  and  reduced  ground  station  complex¬ 
ity.  but  it  requires  that  a  relatively  high  powered  laser  be  flown,  adversely 
impacting  the  space  segment  of  the  system. 

For  this  study,  we  have  assumed  tracking  a  set  of  fixed  point  sources  on  the 
Earth.  For  the  purpose  of  estimating  sensor  performance  and  configuration. 
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active  beacons  are  assumed,  although,  for  the  navigation  simulation,  this 
distinction  is  immaterial. 

This,  then,  is  essentially  a  feasibility  study  for  doing  Autonomous  Navi¬ 
gation  using  a  sensor  of  the  NRL  Rad-Hard  Star  Tracker  class  in  an  Earth- 
based  landmark  configuration. 
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2  AutoNav  System:  Description 

In  this  section,  the  particular  autonomous  navigation  system  studied  is  de¬ 
scribed  in  detail.  The  software  simulation  used  to  test  it  is  described  in 
Sec.  3. 

2.1  Overview 

This  study  examines  the  behavior  of  an  AutoNav  system  comprised  of  a 
single  sensor  very  similar  to  the  NRL  Rad-Hard  star  tracker,  an  Extended 
Kalman  Filter  which  processes  the  sensor  measurements,  and  a  fairly  simple 
Keplerian  orbit  propagator.  A  block  diagram  of  the  system  appears  in  Fig.  1 
on  page  6.  and  the  operation  of  the  system  is  described  below. 

1.  Given  an  initial  estimate,  Uo,  of  the  position  and  velocity  of  the  vehicle, 
the  orbit  propagator  is  used  to  compute  the  position  and  velocity,  u  (t). 
at  subsequent  times.  This  is  a  blind  propagation,  and  is  used  during 
those  intervals  in  which  no  reference  points  fall  within  the  field  of  view 
of  the  sensor. 

2.  When  one  or  more  reference  points  do  fall  within  the  sensor's  ken.  the 
sensor  output,  z,  is  given  to  the  AutoNav  system  for  processing. 

3.  The  AutoNav  system  compares  these  sensor  measurements  to  what  it 
thinks  the  sensor  outputs  should  be.  z.  based  on  its  (probably  erro¬ 
neous)  idea  of  where  the  vehicle  is.  u  ( t ).  The  Kalman  Filter  produces 
the  means  (the  Kalman  gain  matrix,  K)  by  which  these  differences 
in  measurements.  A:,  can  be  converted  into  a  correction,  Au.  of  the 
vehicle's  estimated  position  and  velocity. 

This  process  is  repeated  forever:  fly  blindly  until  measurements  of  known 
reference  points  permit  refinements  of  the  estimated  position  and  velocity 
to  be  made.  Continue  making  these  refinements  until  the  reference  points 
disappear  from  sight,  and  then  go  back  to  flying  blindly. 


U„  CORRECTED  STATE  ESTMATE 

U(t)  PROPAGATED  STATE  ESTIMATE 
AU  OPTIMUM  STATE  CORRECTION 

Figure  1:  AutoNav  System  Block  Diagram 

2.2  Hardware 

2.2.1  Sensor  Description 

The  landmark  sensor  concept  is  based  on  the  design  of  a  slrapdown  star 
tracker  similar  to  ones  currently  in  use.  The  operation  and  performance 
requirements  of  the  two  sensors  are  essentially  identical  and.  as  discussed 
below,  a  similar  set  of  design  tradeoffs  must  be  considered.  A  radiation- 
hardened  version  of  such  a  star  tracker  has  been  built  and  tested  as  part  of 
the  parent  development  program  for  NRL.  That  breadboard  employs  a  rharge 
injection  device  (CID)  detector,  however,  this  sensor  description  also  applies 
to  designs  using  charge  coupled  device  (CCD)  detectors.  The  choice  of  a 
detector  for  this  application  would  depend  on  the  final  system  performance 
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requirements  and  environmental  concerns,  such  as  the  radiation  environment. 
The  state  of  the  art  in  detectors  is  dynamic. 

The  sensor  concept,  shown  in  Fig.  2,  was  developed  as  a  general  purpose 
standard  star  tracker  (REF  1).  A  variation  of  this  design  was  employed  on 
the  Retroreflector  Field  Tracker  (RFT),  which  was  used  to  track  laser  diode 
illuminated  retroreflecting  targets  to  measure  the  motion  of  a  large  flexible 
solar  array  (REF  2).  This  instrument  flew  on  Space  Shuttle  mission  41-D  in 
1984. 


C-RING  GASKET 

CHARGE  INJECTION  DEVICE 
ATE  COOLER  ASSEMLV 


DETECTOR  A  PREAMP  ASSEMBLY 
LOGIC  AREA  ACCESS  COVER 
LOGIC  AREA 

EM  I  FILTERS 

PURER  SUPPLY 
ACCESS  COVER 


OPTICAL  SUPPORT  BASE 
- 19.938  — 


A/N  5160 


Figure  2:  BASCi  Solid  State  Tracker 

t  ypical  performance  and  physical  parameters  for  this  concept  are  shown 
in  fable  2.  A  range  of  parameters  is  shown  to  reflect  that  tradeoffs  in  field  of 
> ■  i  Hi  earn i  ion  lime,  detector  formal,  target  motion,  optic  si/.e.  accuracy, 
etc.  exist  and  t hat  the  sensor  configuration  will  reflect  the  final  requirements. 
1  he  range  of  values  shown  was  used  to  develop  the  landmark  sensor  noise 
models  used  in  the  navigation  filter  simulation. 


i  Field  of  View  <  10  deg  x  10  deg 

I  Accuracy  (temporal)  1  -  10  sec  (ler) 

Target  Rate  0  -  1  deg  /  sec 

Target  Brightness  -8-0  Mv  (equivalent) 

jj  Volume  500  cu  in 

jj  Weight  ,  15  lb  (6.8  kg) 

(  Power _  15  W _ 


Table  1:  Sensor  Capabilities 


The  primary  tradeoffs  are  between  accuracy,  field  of  view,  and  detector 
format.  Constraints  include  target  brightness,  maximum  target  rate,  and 
limitations  on  optics  size. 

The  accuracy  of  the  sensor  defines  the  noise  of  the  landmark  measure¬ 
ment.  Along  with  other  sources,  such  as  attitude  determination  uncertainty, 
it  determines  the  noise  input  to  the  navigation  filter.  This  study  has  con¬ 
sidered  total  temporal  noise  inputs  of  10.  20.  and  30  sec  (1  cr).  In  an  actual 
system  implementation  this  error  would  be  apportioned  to  the  various  sen¬ 
sors  and  subsystems,  but  it  is  apparent  that  a  landmark  sensor  accuracy  on 
the  order  of  that  shown  above  is  adequate. 

Field  of  view  (FOV)  is  a  significant  driver  in  the  sensor  design.  Gener¬ 
ally.  a  small  field  of  view  permits  higher  accuracy  for  a  given  detector  format. 
However,  a  small  field  of  view  limits  the  area  of  the  availability  of  targets,  pro¬ 
viding  less  frequent  updates.  Field  of  view  then  becomes  a  system  trade  that 
is  made  in  the  simulation  of  the  navigation  filter.  This  study  has  considered 
10  and  20  degree  fields  of  view.  Given  typical  available  CTD  detector  formats 
(i.e.  256  to  2048  pixels  square)  a  10  degree  field  of  view  sensor  is  reasonable 
for  the  accuracy  range  considered,  while  a  20  degree  sensor  would  require  use 
of  the  larger  format  detector  and  would  provide  an  accuracy  in  the  higher 
end  of  the  range.  The  Retroreflector  Field  Tracker  used  a  256  x  256  pixel 
CID  with  a  22  x  22  degree  field  of  view.  It  bettered  its  accuracy  specification 
of  19  sec.  Use  of  a  larger  format  detector  would  improve  accuracy. 

The  above  discussion  assumes  a  strapdown.  or  fixed,  landmark  sensor 
with  its  field  of  view  centered  on  nadir.  The  effective  field  of  view  can  be 
increased  by  using  a  one  axis  gimbal  or  a  steering  mirror  to  offset  point  the 
line  of  sight  across  the  satellite  ground  track,  making  more  targets  available. 
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The  current  position  estimate  is  used  to  predict  upcoming  targets  and  the 
sensor  or  line  of  sight  is  offset  to  provide  the  update.  At  low  earth  orbit 
the  horizon  may  be  up  to  70  degrees  from  nadir  and  offset  pointing  would 
provide  significantly  more  updates  from  a  basic  10  degree  sensor.  The  price 
for  this  is  the  increased  complexity  of  the  sensor  and  the  addition  of  pointing 
errors.  Again,  the  tradeoff  between  the  frequency  of  updates  and  accuracy 
is  made  in  the  simulation  of  the  navigation  filter. 

Target  motion  is  an  important  driver  for  the  sensor  design.  Accuracy 
improves  with  the  time  that  the  beacon  signal  is  integrated  on  the  detector. 
Orbital  motion  causes  the  target  image  to  move  on  the  detector,  requiring  the 
track  logic  to  follow  it  and  limiting  the  time  that  a  pixel  integrates  the  signal. 
Maximum  motion  occurs  for  a  target  at  nadir,  and  varies  from  approximately 
1  deg/sec  for  a  450  km  orbit  to  0.0  deg/sec  for  a  geostationary  orbit.  These 
rates  have  been  accommodated  in  star  tracker  designs,  the  typical  tradeoff 
being  in  detector  format:  smaller  formats  produce  a  larger  pixel  angular  field 
of  view  and  allow  a  longer  integration  time,  while  larger  formats  have  greater 
accuracy  due  to  higher  angular  resolution. 

One  difference  between  a  landmark  tracker  and  a  star  tracker  is  that  if 
a  narrow  band  laser  is  used  as  a  beacon  source,  a  narrow  spectral  filter  can 
be  used  to  discriminate  against  the  broad  band  stray  light  background.  This 
technique  was  used  on  the  Retroreflector  Field  Tracker  to  reduce  reflected 
sunlight  from  the  solar  array. 

2.2.2  Beacon  Description 

The  landmark  beacon  configuration  considered  uses  a  near-infrared  (approx¬ 
imately  800  nm  wavelength)  laser  diode  with  an  optical  system  to  shape  and 
project  its  beam.  This  configuration  was  demonstrated  on  the  Retroreflector 
Field  Tracker,  in  which  two  cylindrical  lenses  produced  a  fan-shaped  beam 
which  illuminated  passive  targets  on  the  solar  array.  The  design  requirement 
is  to  provide  a  beam  that  will  produce  a  starlight  level  brightness  over  an 
angular  bear,  width  matched  to  the  sensor  field  of  view. 

Finite  laser  diode  power  levels  provide  a  constraint  on  the  brightness  of 
the  beacon  source.  Current  high  power  laser  diodes  are  available  at  the  Watt 
level.  Table  2  shows  the  brightness  of  a  l-\Yatt  beacon  in  equivalent  star 
magnitude  for  two  beamwidths  and  several  low  earth  orbit  altitudes.  The 
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jj  Width 

Altitude  (km) 

Magnitude 

10  deg 

400 

+4.34 

:i 

'i 

635 

-*-5.34 

;! 

1000 

-*-6.33 

20  deg 

400 

-5.84 

'i 

635 

-6.84 

•i 

1000 

-7.83 

Table  2:  Beacon  Brightness  (Meq) 


beamwidths  (total  angle)  match  the  two  sensor  fields  of  view  used  in  the 
filter  study  and  the  635  km  altitude  is  the  orbit  simulated. 

These  levels  are  within  the  capability  of  the  sensor.  The  use  of  higher 
powered  diodes  or  an  optical  combination  of  several  diodes  or  beacons  would 
increase  the  received  signal.  Each  time  the  power  is  doubled,  the  magnitude 
decreases  (indicating  a  brighter  source)  by  0.75;  for  example,  an  8  Watt 
beacon  with  a  10  degree  beamwidth  would  appear  to  be  a  -3.1  M  object. 
Large  diode  arrays  are  being  developed  as  pumping  sources  for  solid-state 
lasers,  with  powers  up  to  25  Watts  having  been  reported.  If  this  power  were 
used,  the  above  beacon  would  appear  to  be  a  —1.8  M  object. 

Clearly,  a  small  beam  angle  increases  beacon  brightness  and  this  is  one 
factor  that  must  be  considered  in  sizing  the  sensor  field  of  view.  Other 
options  may  be  considered  to  reduce  beacon  beamwidth;  for  instance,  a  fan¬ 
shaped  beam.  A  beam  sized  to  match  the  field  of  view  in  one  dimension 
but  narrower  than  the  field  of  view  in  the  other  would  provide  the  same 
target  availability  for  orbital  tracks  crossing  the  beam  perpendicularly.  The 
target  would  appear,  however,  only  during  a  portion  of  the  pass,  reducing 
the  number  of  times  the  target  could  be  sampled.  Two  such  beams  oriented 
orthogonal  to  each  other  would  provide  availability  for  any  ground  track.  For 
instance,  a  l -Watt,  2deg  xlOdeg  beam  would  appear  as  a  -3.9  M  object, 
instead  of  the  -5.34  M  brightness  shown  above.  An  8-Watt  beacon  would 
appear  as  -1.6  M,  and  a  25- Watt  beacon.  -0.4.  Such  an  improvement  in 
brightness  and  sensor  accuracy  would  have  to  be  traded  against  the  reduced 
number  of  samples  in  t he  filter  simulation. 
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VAX  11/780 

Sun  4/200 

Amdahl  5860 

MIPS  i 

Dhrvstone 

.0  1 

10.9 

16.4 

23.9  ; 

.1 

'  UNPACK 

1.0 

7.9 

— 

25.7  ! 

;  Whetstone 

1.0 

4.0 

— 

16.4  I 

Table  3:  Relative  Performance  Comparisons  of  Processors 


2.2.3  Computer  Hardware 

We  did  not,  in  this  study,  make  an  in-depth  study  of  candidate  flight  com¬ 
puter  systems.  We  merely  recognized  that  the  MIPS  processor,  which  we  are 
targeting  for  other  flight  systems,  is  more  than  equal  to  the  task. 

The  MIPS  is  a  32  bit  RISC  machine  which  can  be  run  at  25  MHz,  but 
is  more  often  operated  at  10  MHz  for  better  matching  to  currently  available 
memory  chips.  It  implements  the  IEEE  754  floating  point  standard,  and 
performs  F.P.  additions  and  multiplications  in  5  clock  cycles,  and  divisions  in 
19  cycles.  It  incorporates  a  5  step  instruction  pipeline,  thereby  attaining  very 
impressive  operation  speeds.  This  is  reflected  in  the  following  comparisons 
to  other  computer  systems  using  several  standard  benchmarks. 

The  figures  in  Table  3  are  for  a  MIPS  running  at  25  MHz,  and  indicate 
that  the  MIPS,  on  rough  average,  has  22  times  the  computing  capabilities  as 
the  VAX.  At  10  MHz.  it  still  outperforms  the  VAX  by  a  factor  of  8.8. 

A  full  compliment  of  program  development  tools  is  also  available  for  the 
MIPS.  The  M120-3  development  station  comprises  a  300  MB  disk  drive.  8 
MB  of  RAM.  and  runs  under  the  full  UNIX  operating  system.  Both  an 
Assembler  and  a  full  optimizing  C  compiler  are  available,  as  well  as  the 
standard,  indispensable.  UNIX  program  development  tools  such  as  grep. 
make,  touch,  and  diff. 

The  CPU  uses  2.5  watts,  and  the  Math  Co-Processor  about  5.0  watts. 
Very  rough  estimates  of  the  power  required  by  a  system  adequate  for  our 
needs  puts  the  power  at  between  12  and  15  watts.  If  power  must  be  limited, 
the  Co-Processor  could  be  eliminated,  and  the  floating  point  computations 
done  in  software  if  the  degradation  in  performance  is  acceptable. 

For  a  fully  contained  MIPS  processing  system  complete  with  boards, 
power  supplies,  and  box.  the  weight  will  be  about  10  pounds. 
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2.3  Algorithms  &  Software 

2.3.1  Propagator 

The  propagator  used  derives  from  the  assumption  that  the  only  force  acting 
on  the  vehicle  comes  from  the  gravitational  field  of  a  homogeneous,  oblate 
Earth.  There  is  no  air  drag,  solar  pressure,  nor  magnetic  effects.  Thus,  we 
have  elliptical  orbits  which  suffer  drift  in  their  ascending  nodes  and  argu¬ 
ments  of  perigee.  The  “state”,  u,  of  the  vehicle  is  taken  to  be  the  standard 
Keplerian  orbital  elements.  That  is, 


a 

e 

i 


UJ 

V 


(1) 


These  terms  are  identified  in  Appendix  H. 

The  differential  equations  of  motion  for  this  choice  of  state  vector  are 
quite  simple. 
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The  propagator  for  this  AutoNav  system  is  merely  a  propagator  for  the 
true  anomaly,  v{t).  Given  a  value  for  t ,  the  true  anomaly  is  computed  in 
several  stages.  First,  compute  the  Mean  Anomaly,  M,  as  follows: 

M  =  —  !(t  -  t0)mod(r)]  (6) 

where  t0  is  the  time  at  epoch,  and  r  is  the  orbit  period.  Incidentally,  the 
quantities  t0  and  r  which  appear  in  the  above  equation  are 


to 


2  tan-1 


e\/l  —  e2  sin  v 
1  ■+■  e  cos  v 


(?) 


and 

r=2*yi.  (8) 

Next,  using  this  value  for  M .  approximate  the  Eccentric  Anomaly,  E,  by 
■‘solving'’  Kepler's  nonlinear  equation 


M  =  E  -  e  sin  E  (9) 

Given  E.  the  True  Anomaly,  v,  is  computed  using  Gauss's  equation. 


v  =  2  tan  ~l 


(10) 


In  the  actual  implementation  of  Eqs.  6,  9,  and  10,  care  must  be  taken  to  keep 
computations  within  the  appropriate  principal  value  bounds.  This  has  been 
done,  but  not  detailed  here  in  the  interest  of  protecting  the  public  welfare. 


2.3.2  Kalman  Filter 

The  basic  equations  of  the  Extended  Kalman  Filter  are  on  display  in  Ap¬ 
pendix  B,  and  the  way  in  which  they  are  used  is  described  there.  Here,  we 
will  merely  present  the  explicit  form  of  certain  of  the  terms  specific  to  this 
problem. 


Linearized  Measurement  Geometry,  H 


It  is  assumed  that  the  output  of  the  sensor(s)  is  preprocessed  to  provide  the 
components  of  a  unit  vector,  given  in  spacecraft  coordinates,  which  points 
from  the  center  of  mass  of  the  vehicle  to  the  target.  If  the  coordinates  of  the 
i-th  reference  point,  in  the  Earth  frame,  £*,  are 

'  *x  ' 

Vi  (11) 

.  z<  . 

then  the  vector,  v^,  from  the  spacecraft  to  this  point  is  given  by 

v,  =  (12) 


<  R3(u;)R1(t)R3(n  -  u )et  —  9o) 


Xi  ' 

a(l  -  e2) 

cos  V 

' 

Vi 

sin  v 

1  —  e  cos  v 

0 

2.  . 

It  is  clear  that  this  is  not  intended  to  be  a  unit  vector.  The  associated  unit 
vector,  g,,  is  simply 

g.  =  (13) 


v,i 


and  is  what  is  produced  bv  the  sensor  system. 

The  matrix  H  (Eq.  47  in  Appendix  B)  is  really  the  concatenation  (stack¬ 
ing  up)  of  the  Jacobian  of  Eq.  13  with  respect  to  u  for  as  many  reference 
points  as  are  visible.  That  is. 


H 


dgi 

3u 


3g> 

au 


(U) 


when  there  are  k  targets  in  view.  In  the  early  stages  of  the  study,  df'^/du 
was  computed  numerically  using  double  sided  finite  differences,  but  this  was 
found  to  somewhat  less  stable  than  desired.  It  is  also  very  computationally 
intensive.  The  Jacobian  is.  therefore,  now  computed  analytically.  This  gives 
the  additional  benefit  of  reducing  the  execution  time  of  the  process.  Actually, 


dg/du  itself  is  not  computed.  Obtaining  these  partials  analytically  would 
be  very  difficult,  even  with  the  aid  of  SMP  .  Instead,  dv/du  is  determined 
analytically  (using  SMP),  and  then  transformed  using  Eq.  77.  repeated  here 
for  convenience. 

<15> 

In  this  equation,  J„  is  the  easily  computed  Jacobian  of  the  non-unit-length 
measurement  vector,  v;  g  is  the  unit  vector  parallel  to  v;  N  is  the  magnitude 
of  v:  I  is  the  unit  matrix;  and  J„  is  the  Jacobian  of  g  . . .  which  is  required 
by  the  Kalman  Filter.  The  analytical  forms  of  v(u)  and  J„  are  given  in 
Appendix  C. 


Linearized  Plant  Geometry,  F 

The  linearized  plant  geometry  matrix,  F  (see  Eq.  46,  Appendix  B)  is  defined 
as 

df 


F  = 


du 


(16) 


where  f(u,t)  is  the  nonlinear  vector  differential  equation  of  motion.  Eq.  2. 
This  matrix  is  also  formed  analytically.  All  its  components  are  zero  except 
for  the  F,,;  listed  below  in  which 
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2Svmbolic  Manipulation  Program;  an  outgrowth  of.  and  improvement  upon.  MIT's 
Macsvma  program  SMP  is  the  creation  of  Steve  Wolfram. 


(21) 


sin  i . 


Fo  =  — r~T 


Fs.i  = 

7(1-5  cos2 i) 

8a 

(22) 

Fs,2  = 

e(l  -  5  cos2  i)  Y 

2(1  -  e=») 

(23) 

FS,3  = 

5  sin  i  cos  i 

2 

(24) 
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State  Covariance  Update,  P*(-) 


The  state  covariance  update  equation.  Eq.  53.  is  used  in  this  study.  It  is  valid 
for  any  gain.  K*.  optimal  or  otherwise,  and  is  reproduced  here  for  reference. 

P*(  +  )  =  {I  -  KfcHfc}Pfc(-){I-  KkHk}r  +  KkRkKTk  (28) 

If  the  optimal  Kalman  gain  is  used  in  Eq.  28.  then  that  equation  reduces  to 

Pfc(-)  -  {I  -  KfcHfe}  Pfc(-)  (29) 

Although  Eq.  29  is  simpler  and  computationally  faster  than  28,  it  is  prone 
to  producing  non-symmetric  Pi,(  +  )  matrices.  This  is  due  solely  to  computa¬ 
tional  roundoff  errors.  Eq.  28,  because  of  its  form,  is  inherently  symmetrical. 
It  is  felt  that  the  freedom  from  worry  about  non-symmetric  updates  is  worth 
the  additional  computation  load  of  using  Eq.  28.  Reducing  computation 
speed  is  not  a  high  priority,  right  now,  anyway.  When  this  algorithm  is  com¬ 
mitted  to  flight  code,  then  it  will  be  worth  examining  this  matter  to  see  if 
Eq.  29  can  be  safely  used. 
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2.3.3  Computer  Languages 

The  software  for  this  entire  project  was  written  in  the  C  language  using  the 
Microsoft  C  Compiler,  Version  5.1.  This  language  was  chosen,  primarily,  for 
the  following  reasons: 

1.  It  is  ideally  suited  for  scientific  applications,  and  eclipses  FORTRAN 
in  these  areas. 

2.  A  vast  store  of  support  software  is  available  within  the  C  community 
which  greatly  aids  and  simplifies  the  software  development  process.  No 
other  language  enjoys  this  support. 

3.  C  was  designed,  at  its  inception,  to  be  both  an  easily  portable  language, 
and  to  produce  efficient  executable  code.  Because  of  this,  the  code 
written  for  this  project  can  be  very  easily  moved  onto  a  flight  processor. 
BASG  has  experience  in  this  area.  We  extracted  the  computational 
kernel  of  CLOP3  for  use  as  the  real-time  orbit  propagator  in  one  of 
our  satellites4.  The  language  is  equally  well  suited  for  scientific  and 
real-time  applications. 

2.3.4  Reference  Points 

The  coordinates  of  the  reference  points  (targets)  used,  in  this  study,  for  the 
Kalman  update  of  the  estimated  state  are  contained  in  a  list  of  89  Air  Force 
bases  located  in  the  political  U.S.A.  The  east  longitude,  north  latitude,  and 
height  above  sea  level  for  each  base  is  specified.  The  file  of  these  coordinates 
is  given  in  Appendix  G. 

It  is  assumed  that  each  of  these  bases  houses  a  laser  which  can  be  seen 
by  the  sensor.  See  Sec.  2.2.2  for  a  fuller  treatment  of  this  subject. 

REFERENCES: 

3A  general  purpose,  numerically  integrating,  orbit  propagator  developed  at  BASG. 
written  by  Charlie  Ros5e,  comprising  high-order  gravitational  potentials,  a  variety  of  air- 
drag  models,  and  designed  to  run  on  personal  computers. 

'The  Remote  Mirror  Experiment,  R.V1E. 
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3  AutoNav  System:  Evaluation 

In  this  section,  we  describe  how  the  AutoNav  system  was  evaluated,  and  what 
simplifications  were  instituted  in  order  to  finish  the  study  in  finite  time. 

3.1  Overview 

The  AutoNav  system,  and  the  vehicle  within  which  it  lives,  necessarily  inter¬ 
act  with  the  universe.  The  satellite  is  acted  upon  by  the  gravitational  fields 
of  the  Earth.  Moon,  Sun,  planets.  It  is  buffeted  about  by  solar  pressure  (both 
direct  and  reflected  by  the  Earth),  air  drag,  and  magnetic  fields  for  starters. 
The  sensors  produce  noisy  measurements,  through  a  nonlinear  transfer  pro¬ 
cess.  of  reference  points  viewed  through  a  shimmering,  refractive  medium. 
The  final,  and  ultimate,  evaluation  of  any  AutoNav  system's  operation  must 
take  place  in  this  environment. 

Such  real  world  tests,  obviously  and  necessarily,  are  prohibitively  expen¬ 
sive.  They  also  do  not  permit  experimentation  with  the  algorithms.  What 
is  done,  at  this  stage  of  the  investigation,  is  to  simulate  the  universe  with  as 
much  fidelity  as  one  can  afford.  The  economics  of  the  situation  involve  time, 
money,  computing  facilities,  and  the  desired  accuracies  of  the  the  results. 
A  block  diagram  of  the  simulation  used  to  evaluate  the  AutoNav  system  is 
shown  in  Fig.  3. 

In  a  comprehensive  evaluation,  the  following  quantites  would  be  simu¬ 
lated.  Such  a  complete  simulation  was  not  warranted  at  this  stage  of  our 
investigations,  and  the  areas  of  departure  are  discussed  in  the  following  sub¬ 
sections. 

•  Forces  on  the  vehicle.  As  mentioned  earlier,  they  include  those  arising 
from  the  gravitational  fields  of  many  sources:  air  drag,  magnetic  forces 
if  the  vehicle  is  charged,  and  solar  pressures  both  directed  and  reflected. 
The  simulation  of  these  forces  is  embodied  in  what  we  term  the  “Truth 
Propagator”.  This  propagator  must  represent  the  actual  motion  of  the 
vehicle  as  closely  as  possible,  and  is  in  contrast  to  the  simpler  "AutoNav 
Propagator”  which  forms  part  of  the  AutoNav  system. 

•  The  sensors.  Measurement  noises,  nonlinearities.  and  mounting  mis¬ 
alignments  must  be  modeled. 
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Figure  3:  Simulation  for  Auto.Nav  System  Evaluation 

•  Tlie  reference  points.  Uncertainties  in  their  surveyed  positions  may 
he  modeled,  as  well  as  uncertainties  in  their  observed  positions  due  to 
atmospheric  effects.  These  uncertainties  comprise  both  random  fluctu¬ 
ations.  and  slowly  drifting  biases. 

•  The  vehicle  attitude.  This  is  important  only  for  t  hose  AutoXav  systems 
which  make  angular  measurements  upon  the  reference  points.  Ranging 
schemes  are  usually  not  plagued  with  such  complications. 

3.2  Landmarks 

It  was  deemed  appropriate,  and  realistic,  to  choose  the  Air  Force  bases  lo¬ 
cated  in  the  political  United  States  as  the  locations  of  the  reference  points 


|  Orbit 
!  Parameter 

Value  i 

i 

j  a 

7000.000  km 

■1  e 

0.010  i 

J  i 

54°40'  j 

ft 

359.984°  ■! 

UJ 

0.011°  i 

V 

56.413°  i 

Table  4:  Definition  of  Test  Flight  Orbit 

which  the  star  tracker  would  observe.  Bases  located  in  other  parts  of  the 
world,  especially  Europe,  would  be  the  first  to  become  unavailable  at  the 
precise  time  when  the  functioning  of  the  AutoNav  system  became  critical. 
The  latitude,  longitude,  and  height  above  sea-level  of  89  such  bases  were 
placed  in  a  data  file,  and  accessed  by  the  simulation  program.  A  listing  of 
this  data  file  is  included  in  Appendix  G.  The  locations  of  the  bases  are  shown 
in  Fig.  4  on  page  58. 

Recall  that  the  actual  targets  being  observed  are  laser  beacons  located 
at  these  Air  Force  bases.  To  simplify  the  building  of  the  simulator,  it  was 
assumed  that  each  of  these  beacons  radiated  uniformly  over  a  locally  vertical 
cone  having  the  same  field  of  regard  as  the  tracker.  Thus,  so  long  as  the 
beacon's  location  was  in  the  field  of  view  of  the  tracker,  the  tracker  responded. 
This  is  an  over-simplification,  and  will  be  corrected  in  subsequent  studies. 

Uncertainties  in  the  surveyed  positions  of  the  landmarks  will  be  manifest 
as  small  angular  offset  biases  on  the  order  of  10  sec,  and  must  be  consid¬ 
ered.  However,  they  were  ignored  in  this  study  because  of  the  gross  errors 
introduced  by  the  AutoNav  propagator.  They  will  not  be  omitted  in  the 
future. 


3.3  Orbit 

An  orbit  with  the  following  initial  conditions,  given  as  Keplerian  elements, 
was  chosen  as  being  typical  of  the  flights  that  could  be  expected. 

The  right  ascension  of  the  ascending  node,  ft.  was  chosen  to  place  the 
ground  track  of  the  orbit  over  the  most  A.F.  bases  along  the  east  coast.  This 
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gave  us  the  most  measurements  possible:  an  important  consideration  given 
the  difficulties  we  experienced  with  our  choice  of  Model  Propagator5.  The 
ground  track  for  this  orbit  is  shown  in  Fig.  4.  The  locations  of  the  A. F. bases 
are  denoted  by  -  signs. 

3.4  Tracker  Configuration 

The  AutoNav  system  included  a  single  tracker  having  the  characteristics 
described  in  Sec.  2.2.1.  It  was  aligned  to  nadir  to  within  the  accuracy  capa¬ 
bilities  of  an  assumed  attitude  control  system.  This  ACS  was  not  modeled; 
its  errors  were  assumed  to  be  zero-mean  with  a  Gaussian  distribution,  and 
w'ere  simply  RSSed  with  the  intrinsic  beacon  tracker  noise.  The  same  is  true 
for  other  uncertainties  in  the  measurement  systems,  such  as  ‘Target  twinkle'1 
arising  from  atmospheric  instabilities. 

Constant,  or  very  slowly  varying,  biases  were  ignored  in  this  study.  These 
included  sensor  mounting  offsets,  and  the  effects  of  atomospheric  refraction. 
In  these  preliminary  investigations,  it  was  assumed  that  refraction  can  be 
predicted,  and  appropriate  compensations  made. 

The  field  of  view  of  the  tracker  could  be  specified  as  a  run  time  parameter. 
Flights  with  different  FOVs  were  made  so  it  could  be  determined  if  a  gimbaled 
mirror  scanner  would  have  to  be  added  to  the  sensor  system  in  order  to  see 
enough  targets. 

The  tracker  noise  was  also  considered  to  be  Gaussian  with  zero-mean. 
The  noise.  rp.  was  simulated  by  using  the  following  simple  approximation  to 
a  Gaussian  generator: 

12 

rp  =  fir  *  <Tr  -  °’5)  (3°) 

fc= i 

where  <rr  and  /i,  are  the  desired  standard  deviation  and  mean  of  the  pseudo¬ 
random  sequence.  { rp },  and  u  is  a  random  variable  which  is  uniformly  dis¬ 
tributed  on  the  interval  0.  1  .  Granted:  this  is  a  kindergarten  variety  of 
Gaussian  generator,  but  it  was  deemed  adequate  because  our  primary  con¬ 
cerns  lay  in  the  gross  behavior  characteristics  of  the  system. 

5See  Sec.  F  where  this  problem  is  discussed  in  detail 


3.5  Initial  Covariance  Matrices 


The  Kalman  Filter  incorporates  three  covariance  matrices,  the  values  of 
which  affect  the  operation  of  the  system.  See  Appendix  B  for  a  fuller  de¬ 
scription  of  these  matrices. 

3.5.1  State:  P 

This  is  a  time  varying  matrix  whose  value  is  propagated  by  solving  the  dif¬ 
ferential  equation.  Eq.  50,  and  corrected  occasionally  by  Eq.  53.  Eq.  50  must 
be  given  an  initial  value  for  P  at  the  beginning  of  the  simulation,  and,  in  lieu 
of  good  reasons  for  choosing  otherwise,  the  following  was  used: 
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3.5.2  Plant:  Q 


Often,  the  Q  matrix  is  set  to  0  for  simplicity  of  operation,  and  because  the 
analyst  has  no  theoretical  basis  for  choosing  anything  different.  The  result 
of  this  is  that  the  gain  matrix.  K,  goes  to  zero  after  a  while,  and  the  Kalman 
Filter  becomes  "conceited”  by  refusing  to  produce  corrections  to  the  estimate 
of  the  state. 

We  chose  Q  to  be  the  following  based  on  our  desire  to  force  the  derivative 
of  P  to  have  a  certain  non-zero  value  when  P  was  zero.  This  is  what  we 
used: 


10"'  0  U  0  0  0 

0  10"7  U  0  0  0 

0  0  I0"7  000 

000  10"9  0  0 

o  o  o  a  nr 9  o 

0  0  0  0  0  I0"9 


(32) 
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3.5.3  Measurement:  R 


The  dimension  of  R  varies  with  the  number  of  measurements  being  made. 
Therefore,  it  could  not  be  specified  as  a  particular  thing  at  run  time.  We  did 
define  R,  at  each  iteration,  to  be  a  constant  diagonal  matrix  whose  diagonal 
elements  were  the  variance,  <t2,  of  the  tracker  noise. 


3.6  Truth  Propagator 

The  Truth  Propagator  was  initially  intended  to  be  only  the  output  of  CLOP 
using  a  set  of  spherical  harmonic  zonal,  sectoral,  and  tesseral  coefficients 
of  fairly  high  order.  Although  this  was  done,  we  found  it  necessary  to  use 
another,  simpler  propagator  to  verify  the  correct  operation  of  the  AutoNav 
system.  This  is  described  below. 

Although  CLOP  has  the  capability  of  including  air-drag  forces,  we  flew  all 
our  simulations  in  vacuum.  This  is  because  the  AutoNav  propagator  errors 
overshadowed  all  others,  and  air-drag  effects  would  have  been  unnoticable. 

3.6.1  High-Order  Harmonic 

The  coefficients  used  to  represent  Absolute  Universal  Truth  were  obtained 
from 
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Although  coefficients  through  36th  order  are  given  there,  we  used  only  those 
through  order  16.  W'e  estimated  that  the  effects  of  higher  order  terms  would 
be  miniscule  for  the  flight  durations  we  were  considering. 
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3.6.2  Simple  Truth  Propagator 


Attention!  This  subsection  describes  a  simple  test  we  made  to  verify  the 
correct  operation  of  the  simulation.  It  served  no  other  function  than  that! 

Once  we  realized  that  an  actual  orbit  differed  substantially  from  one  rep¬ 
resented  by  Keplerian  elements  in  which  f!  and  ui  drifted  at  constant  rates 
given  by  Eq.  3  and  4,  we  temporarily  replaced  the  accurate,  CLOP  propaga¬ 
tor  with  the  Keplerian  AutoNav  propagator.  This  permitted  us  to  evaluate 
system  performance  in  the  absence  of  divergence  between  propagators,  to 
validate  the  Kalman  Filter,  and  to  gain  basic  insight  into  the  operation  of 
the  system  as  a  whole. 

The  filter  brought  the  initial  error  from  a  kilometre  down  to  a  millimetre 
very  quickly,  thereby  giving  us  confidence  in  the  correctness  of  the  system  as 
a  whole. 

3.7  System  Performance 

We  evaluated  the  AutoNav  system  performance  by  making  a  series  of  sim¬ 
ulated  flights  over  the  eastern  United  States.  This  ensured  that  the  sensor 
would  see  enough  targets  to  keep  the  Truth  and  AutoNav  propagators  from 
diverging  too  much.  Plots  of  the  position  errors  for  the  various  flights  ap¬ 
pear  figures  6  through  17.  and  the  conditions  of  the  flights  are  summarized 
in  the  following  table.  The  reader  is  cautioned  that  the  vertical  axis  scales 
on  these  plots  are  almost  all  different.  We  accepted  the  automatic  scaling  of 
the  graphing  package  to  show  the  maximum  detail  in  the  plots. 

In  the  Truth  Model  column  of  the  following  tabl.e.  the  term  "Keplerian'' 
means  that  the  Approximate  Truth  propagator  described  in  Sec.  3.6.2  was 
used.  The  term  ‘‘16th''  denotes  the  16th  order  spherical  harmonic  model 
based  on  the  LAGEOS  gravitational  potential  model. 

The  position  error  plots  for  the  various  flights  are  probably  self-explanatory, 
but  a  few  observations  can  be  made.  The  reader  will  likely  notice  that,  in 
many  cases,  even  though  the  position  error  has  been  reduced  to  very  low  lev¬ 
els.  the  position  error  begins  to  increase  substantially  when  measurements 
are  no  longer  available.  This  is  because  the  velocity  error  results  in  an  es¬ 
timated  orbit  which  has  a  slightly  different  inclination  than  the  real  orbit. 
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FOV  l<r  Noise 

Truth  Model  Init.  Error  1 

10°  0  sec 

Keplerian  100  m  ! 

20°  0  sec 

Keplerian  1  100  m  i 

Keplerian 


Keplerian 


16th 


16th 


100  m 


2000  m 


0  m 


0  m 


0  m 


|  lb  I  20°  0  sec 


lc  20°  0  sec 


2  10°  0  sec 


■;!  2a  |  20°  0  sec 


j  10°  10  sec  !  16th  |  0  m  i 


'  3a  |  20°  10  sec 


•\  3b  I  10°  20  sec 


2 


;|  3d  10°  30  sec 


j!  3e  20°  30  sec 


Table  5:  Synopsis  of  Flight  Parameters 

Hence,  the  estimated  position  begins  to  draw  away  from  the  real  position. 

Small  u  Error  Fig.  5  shows  how  the  distance  between  two  vehicles  changes 
when  their  initial  positions  and  velocities  are  slightly  different.  For  the 
case  shown,  the  position  error  was 


16th 

- o^TT 

16th 

0  m  j 

16th 

0  m  i 

16th 

0  m  | 

16th 

0  m  j 

Ap  = 


metres 


Av  = 


metres  sec  1 


Note  the  axes’  units.  The  vertical  axis  is  in  kilometres,  while  the  hor¬ 
izontal  axis  is  in  thousands  of  seconds.  Thus,  after  some  three  orbits, 
this  very  small  initial  error  grown  to  about  6  kilometres. 

Kepler  Representation  Error  The  plot  shown  in  Fig.  18  on  page  65  shows 
how  the  position  of  a.  vehicle,  as  predicted  by  the  Keplerian  formula¬ 
tion  with  drifting  f!  and  u>.  departs  from  an  actual  orbit  with  which 
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it  initially  matches.  In  this  case,  the  ‘‘actual  orbit”  is  derived  bv  us¬ 
ing  a  sixteenth  order  geopotential  model  to  compute  the  gravitational 
forces,  an  all  other  forces  are  ignored.  For  a  near-circular  orbit  having 
a  54.666°  inclination,  and  100  min  period,  the  Keplerian  prediction  is 
in  error  by  some  16  km  after  a  single  orbit.  This  indicates  that  it  is  a 
poor  candidate  for  an  AutoNav  propagator. 

Flight  1  The  AutoNav  system  was  given  a  100  m  initial  error.  Forty  three 
measurements  on  several  targets  were  made  during  the  flight,  and  it 
is  seen  that  the  error  was  continually  reduced  until  no  more  targets 
could  be  seen.  Then,  since  the  Truth  and  AutoNav  orbits  were  slightly 
displaced  despite  their  identical  analytical  form ,  the  true  and  estimated 
positions  began  to  diverge.  (Fig.  6) 

Flight  la  This  is  the  same  as  Flight  1,  except  that  the  sensor's  FOV  is 
20  degrees.  Now  264  measurements  were  made  on  about  ten  different 
beacons,  and  the  initial  error  was  reduced  substantially.  The  very  slight 
residual  error  remaining  after  about  275  seconds  was  still  enough  to 
cause  ultimate  divergence.  (Fig.  7) 

Flight  lb  This  is  the  same  as  Flight  la,  but  shows  behavior  over  a  full  orbit. 

rather  than  just  the  time  during  which  measurements  are  available.  (Fig,  8) 

Flight  lc  Again,  essentially  the  same  as  Flight  1,  but  with  an  initial  position 
error  of  2  kilometres.  The  purpose  of  this  flight  was  to  show  that  the 
Kalman  Filter  could,  indeed,  correct  the  AutoNav  propagator  after  the 
system  had  flown  with  measurements  for  a  full  orbit.  (Fig.  9) 

Flight  2  Here,  the  16th  order  LAGEOS  gravitational  potential  was  used 
for  the  truth  model.  The  large  difference  between  it  and  the  AutoNav 
Keplerian  propagator  is  evident  in  the  very  large  error  growth  rates: 
roughly  3.5  metres/ sec  when  no  measurements  are  available.  (Fig.  10) 

Flight  2a  The  same  as  Flight  2.  but  with  a  20  degree  sensor  FOV'.  Again: 
large  error  growth  rates  are  seen  in  the  absence  of  measurements  but, 
because  of  the  larger  aspect  angle.  265  measurements  were  made  which 
kept  the  maximum  error  lower.  (Fig.  11) 


The  remaining  flights  were  made  with  the  16th  order  LAGEOS  model,  and 
the  matrix  of  FOVs  and  sensor  noises  shown  in  the  table.  They  were  intended 
to  give  us  insight  into  the  sensor  and  attitude  control  system  uncertainties 
that  could  be  tolerated  for  acceptable  navigation  performance.  Unfortu¬ 
nately,  the  inadequacy  of  the  AutoNav  propagator  caused  estimation  errors 
which  overshadowed  the  effects  of  the  sensor  noise. 


4  Conclusions 

The  question  of  whether  a  beacon  tracker,  similar  to  the  NRL  Rad-Hard 
Star  Tracker,  can  be  successfully  used  in  a  Landmark  Tracking  AutoNav 
application  rests  squarely  on  the  accuracies  required,  the  number  of  beacons 
available,  and  their  optical  characteristics  as  well  as  distribution  over  the 
surface  of  the  Earth.  Basically,  we  have  found  that  a  star  tracker  can  be  used 
in  a  landmark  approach,  but  only  within  some  rather  severe  constraints: 

•  Beacons  can  be  distributed  more-or-less  uniformly  over  the  Earth.  Re¬ 
striction  to  the  United  States  does  not  seem  feasible." 

•  Position  errors  of  several  kilometres  can  be  tolerated  between  updates.' 

4  A  very  sophisticated  orbit  propagator  is  built  into  the  system. 

There  are  non-trivial  problems  associated  with  this  approach  to  the  problem. 
Most  of  them  stem  from  the  fact  that  the  tracker  would  be  used  in  a  way  for 
which  it  was  not  designed  nor  intended. 

•  The  beacons  "twinkle”/ 

•  The  light  from  the  beacons  is  refracted  by  the  atmosphere  by  varying 
amounts  which  can  be  only  partially  compensated  for." 

•  The  tracker's  field  of  view  is  too  small:  it  greatly  reduces  the  number 
of  targets  that  can  be  seen,  as  well  as  the  durations  of  observations. 
This  results  in  fewer  updates,  and  necessitates  a  very  sophisticated 
orbit  propagator."  If  the  FOV  is  increased  optically,  the  measurement 
errors  increase  in  proportion.  If  it  is  expanded  by  a  gimbaled  "steering 
mirror”,  the  mechanical  complexity  is  increased  greatly,  as  is  the  system 
weight.  Accuracies  also  suffer.  When  the  FOV  is  enlarged  by  any 
means,  there  is  an  attendant  increase  in  the  complexity  of  the  beacons 
themselves.  Their  field  of  radiation  must  be  enlarged  accordingly. 

•  The  already  too  few  beacons  can  be  obscured  bv  cloud  cover.  This  is. 
perhaps,  the  most  egregious  problem  with  this  approach. 


4.1  Why  The  Keplerian  AutoNav  Propagator? 

The  question  naturally  arises  as  to  why  we  chose  to  use  the  Keplerian  for¬ 
mulation,  with  fi  and  ur  drifting  at  the  rates  predicted  by  the  First  Order 
Theory,  for  the  AutoNav  propagator.  It  was  done  for  two  reasons. 

Simplicity  A  guiding  principle  in  designs  of  this  kind  is  to  make  the  problem 
as  simple  as  possible.  By  choosing  the  Keplerian  representation,  rather 
than  the  Cartesian,  we  hoped  to  reduce  the  computational  load  on 
the  flight  computer,  speed  up  the  operation,  and  reduce  the  size  and 
complexity  of  the  hardware. 

Ignorance  We  did  not  realize,  in  the  beginning,  that  the  actual  orbital 
motion  would  depart  so  drastically  from  the  First  Order  Theory  as  it 
did.  W'e  knew  that  there  would  be  some  difference  between  reality  and 
the  model,  but  did  not  suspect  that  it  would  be  this  great.  If  we  had, 
we  would  have  gone  to  the  Cartesian  representation  (with  its  ability  to 
easily  handle  complex  forces  from  various  sources)  at  the  outset,  and 
accepted  the  cost  of  computational  complexity. 

4.2  Suggestions  for  Further  Study 

With  some  enhancement  of  already  existing  simulation  software,  we  will  be 
able  to  answer,  with  confidence,  the  questions  of  accuracies  and  hardware 
requirements  of  an  optical,  landmark  tracking  AutoNav  system.  Such  en¬ 
hancements  would  include  the  following: 

•  Replace  t he  simplistic  Keplerian  propagator  in  the  AutoNav  system 
with  one  which  solves  the  differential  equations  of  motion  numerically. 
This  involves  only  lifting  the  computational  kernel  from  CLOP,  which 
already  exists.  This  propagator  would  be  formulated  in  terms  of  Carte¬ 
sian  coordinates,  not  the  Keplerian  orbital  elements.  (See  Sec.  F) 

•  Casting  the  AutoNav  propagator  in  terms  of  Cartesian  coordinates  will 
necessitate  reformulating  the  linearized  plant  and  measurement  matri¬ 
ces.  F  and  H,  for  the  Kalman  Filler.  This  is  not  difficult,  and  will 
result  in  simpler  computations.  We  also  suspect  that  somewhat  lower 
state  estimation  errors  will  result  from  such  a  formulation,  but  it  is 
only  that:  a  suspicion. 
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•  As  a  means  of  refining  the  Truth  Model,  add  a  sophisticated  model 
of  air  density  to  CLOP.  This  will  probably  be  one  of  the  Jacchia- 
Roberts  models,  with  which  we  already  have  some  familiarity.  These 
are  dynamic  models  which  permit  the  density  to  vary  with  the  relative 
position  of  the  Sun. 

•  Add  the  position  of  the  sun  to  the  simulator  truth  model.  This  will 
permit  day-night  discriminations  to  be  made  when  deciding  if  a  given 
target  is  visible.  The  software  module  to  do  this  already  exists. 

•  Include  bias  errors  such  as  sensor  mounting  position,  sensor  signal  off¬ 
sets,  reference  point  location,  and  atmospheric  refraction  in  the  para¬ 
metric  studies. 

•  Perform  a  comprehensive  sensitivity  analysis  of  the  system  which  will 
tell  us  the  effects  of:  simplifying  the  AutoNav  gravitational  potential 
model,  air  drag,  bias  errors  in  various  parts  of  the  system,  etc. 

•  Recast  the  entire  AutoNav  system  in  terms  of  using  range,  rather  than 
angular ,  measurements  to  beacons  located  either  on  the  ground  or  on 
cooperative  satellites  in  (probably)  synchronous  Earth  orbit.  We  feel, 
for  a  variety  of  reasons,  that  such  an  approach  would  provide  a  sim¬ 
pler,  more  accurate,  and  more  reliable  solution  to  the  DOD  AutoNav 
problem. 


31 


This  Page  Left  Blank  On  Purpose 


32 


Append  ices 
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A  Notation 


The  symbol  E  is  used  to  represent  an  orthogonal,  right-handed,  coordinate 
system  in  3-space.  Superscripts  distinguish  one  system  from  another,  while 
subscripts  denote  coordinate  axes.  Thus,  Eq  specifies  the  q  coordinate  sys¬ 
tem,  while  E%  denotes  the  2-axis  of  coordinate  system  q.  The  subscripts  can 
only  have  values  1,  2,  and  3. 

Vectors  are  denoted  by  boldface  lower  case  letters,  while  matrices  are 
boldface  upper  case.  Thus,  the  product  of  vector  v  by  matrix  M  to  produce 
vector  z  is  written  as 

z  =  Mv  (35) 

The  term  R*(0)  represents  an  orthogonal  transformation  (rotation)  of  a 
coordinate  frame  about  its  k  axis  through  an  angle  0.  Such  rotations  are 
right-handed.  Thus, 

'1  0  O' 

Ri(0)  =  0  cos(0)  sin(0)  (36) 

0  -  sin(0)  cos(0) 

cos(0)  0  -  sin(0) 

0  1  0  (37) 

sin(0)  0  cos(0) 

cos(0)  sin(0)  0 

-  sin(0)  cos(0)  0  (38) 

0  0  1. 

So.  if  coordinate  system  Eb  is  produced  by  rotating  coordinate  system  Ea 
about  Ea' s  2-axis  through  an  angle  3.  this  transformation  is  expressed  as 

Eh  --  R 2(J)E°  (39) 

Next,  suppose  we  have  a  vector,  r.  whose  components  in  Ea  are  given  as  r°. 
Its  components  in  Eh  may  then  be  computed  as 

r-  -  R2(J)r°  (40) 

A  sequence  of  rotations  is  expressed  as  an  ordered  sequence  of  products 
of  these  rotation  matrices.  For  example,  a  rotation  about  the  1-axis  through 


R2(0)  - 

R3(0)  = 
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an  angle  followed  by  a  rotation  about  the  new  3-axis  through  and  angle 
02  is  written  as 

T  =  R3(02)R1(01)  (41) 

Remember  that  the  order  of  these  terms  is  crucial! 

An  obvious  identity  obtains  from  the  fact  that  two  sequential  rotations 
about  the  same  axis  is  the  same  as  a  single  rotation  about  that  axis  through 
the  sum  of  the  individual  angles.  Thus. 

R*(a)R*(J)  =  R*(a  -  3)  (42) 

Since  rotation  matrices  have  the  special  property  that 

(Rk(oO)"1  =  Rfc(-a)  (43) 

the  inverse  of  a  sequence  of  rotations  can  be  expressed  simply  bv  reversing 
the  order  of  mulitplication,  and  changing  the  signs  of  their  angles.  That  is. 

(R7(a)Rfc(  •?))-'  =  R*(-tf)Rj(-a)  (44) 

and  so  forth. 
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B  Extended  Kalman  Filters 


This  section  is  merely  a  presentation,  and  a  minimal  elucidation,  of  the 
equations  of  the  Extended  Kalman  Filter.  No  attempt  to  derive  the  equations 
is  made.  If  the  reader  is  interested  in  detail  at  that  level,  there  are  a  multitude 
of  references  where  s/he  may  look.7 

The  following  terms  appear  in  the  Kalman  equations: 


Symbol 

Description 

Dimension 

n 

N  umber  of  state  vector  components 

m 

Number  of  measurements 

u 

State  vector 

n 

z 

Measurements  vector 

m 

p 

State  covariance 

n  x  n 

R 

Measurements  covariance 

m  x  m 

Q 

Plant  covariance 

Tl  X  71 

K 

Kalman  gain 

n  x  m 

F 

Linearized  plant  geometry 

n  x  n 

H 

Linearized  measurement  geometry 

m  x  n 

The  following  compact  notation  is  also  used: 


q{~)  =  q{tk  x  St)  (45) 

where  St  is  an  infinitesimal  time  interval,  and  the  subscript  fc  denotes  a 
variable's  value  at  the  A:-th  sample  time. 

Let  h(u.f)  be  the  (probably)  nonlinear  measurement  prediction,  and  let 
f(u./)  be  the  (almost  assuredly)  nonlinear  differentia!  equations  of  motion. 
Furthermore,  define  the  Jacobians  of  f  and  h  to  be 


di 

du 


(46) 


H 


dh 

cfu 


(17) 


See  for  instance.  "Applied  Optimal  Optimization".  Arthur  Gelb,  Editor.  MIT  Press. 
1071  Chapter  fi.  page  180.  treats  nonlinear  estimation 
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That  is, 


(48) 


duk 


and  similarly  for  H. 

During  periods  when  no  measurements  are  available,  the  state  vector  and 
the  state  covariance  matrix  are  propagated  forward  in  time  by  solving  the 
following  differential  equations: 


du 

It 


=  f(u.t) 


(49) 


^  =  FP  -  PFr  -t-  Q(u,t) 


These  two  differential  equations  are  solved8  in  the  interval 


(50) 


tk- 1  <  t  <  tk 


(51) 


At  time  4,  when  a  measurement  set  is  available,  those  measurements  are 
processed  using  the  following  Kalman  equations: 

K*  =  P*( -)H[{H*P*(-)H[ -R*}"1  1 52) 

P*( -)  =  {I-  KfcHfc}Pfc(-){I  -  KkUk}T  -  KkRkKTk  (53) 

Ufc(~)  =  Ufc(-)  -  Kfc{xfc  -  h(ufc(-).t)}  (54) 

Equations  53  and  54  are  then  used  as  initial  conditions  for  49  and  50 

over  the  next  dry  spell  of  no  measurements. 


'‘Solved.  in  all  likelihood,  bv  using  numerical  algorithms  such  as  Runge-Kutta  or 
Bulirsch-Stor 
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C  Analytical  Expressions  for  v  and  dw / d\i 

This  appendix  presents  the  analytical  expressions  for  v(u)  and  d\ id u.  They 
have  been  lifted,  bodily,  from  the  SMP  log  files  in  which  they  were  derived. 
They  are,  therefore,  functionally  correct,  but  not  pretty. 

The  following  are  a  set  of  replacement  rules  (variable  name  substitutions) 
that  are  pertinent  to  the  rest  of  the  appendix.  It  should  be  noted  that 


Omega 

— ♦ 

n 

omega 

— ♦ 

aJ 

se 

— ♦ 

t 

- > 

t 

i 

— ♦ 

i 

nu 

— ♦ 

1/ 

Co  si] 

— ♦ 

cos() 

Sin;] 

— ♦ 

sin() 

v'kj 

— ♦ 

Vfc 

jvij.kj 

— * 

(Cos [Omega  -  t*we 
(Sin [Omega  -  t*we 
(Cos [mi  +  omega] 
(Sintnu  +  omega] 
(Cos  [i] 

(Sin[i] 

(Cos [nu] 

(Sin[nu] 

(Cos [omega] 

(Sin [omega] 


thetaO] 

thetaO] 


->  cl) 
->  si) 
->  c2) 
->  s2) 
->  c3) 
->  s3) 
->  c4) 
->  s4) 
->  c5) 
->  s 5 ) 


(sl*s2  -  cl*c2*c3 
(cl*s2  +  c2*c3*sl 
(cl*c2  -  c3*sl*s2 
(c2*sl  +  cl*c3*s2 


->  terml) 
->  term2) 
->  term3) 
->  term4) 


The  vector,  v,  from  the  vehicle  to  the  reference  point  at  x,  y,  z'T  is  given 
by 

v[l]  :  -(terml*y  +  term2*x  -  (c2*s3*z)) 
v[2]  :  -(c3*z  -  (s3*(cl*y  -  (sl*x)))) 

v[3]  :  (a*(l  -  a'2))/(l  +  c4*e)  -  (tenn3*x)  -  (term4*y)  - 
(s2*s3*z) 

As  noted  earlier,  this  is  not  a  unit  vector. 

The  components  of  the  Jacobian,  dv/du,  are 
jvCl.l]  :  0 
jv[l  ,2]  :  0 

jv[l,3]  :  c2*(c3*z  +  s3*(-(cl*y)  +  sl*x)) 

jv[l,4]  :  terml*x  -  (tenn2*y) 

j v [1 , 5]  :  -(term3*x  +  terni4*y  +  s2«s3*z) 

jv[l,6]  :  -(term3*x  +  tenn4*y  +  s3*z*(c4*s5  +  c5*s4)) 

j  v  [2 , 1]  :  0 

jv[2 ,2]  :  0 

jv[2,3]  :  c3*(cl*y  -  (sl*x))  +  s3*z 
jv [2 ,4]  :  s3*(-(cl*x)  -  (sl*y)) 
j  v  [2 , 5]  :  0 


jv[2,6]  :  0 

jv[3,l]  :  (1  -  a*2)/(l  +  c4*e) 

jv [3,2]  :  - ( (a*(2e  +  c4*(l  +  e‘2)))/(l  +  c4*a)‘2) 

jv[3,3]  :  -(s2*(c3*z  +  s3*(-(cl*y)  +  sl*x))) 

jv[3,4]  :  -(term3*y)  +  term4*x 

jv[3,5]  :  terml*y  +  term2*x  -  (c2*s3*z) 

jv[3,6]  :  -( (a*a*s4*(-l  +  e‘2))/(l  +  c4*a)*2)  +  tennl*y  + 
term2*x  -  (c2*s3*z) 


Hi 


D  Transformation  of  Jacobians 


Here  is  the  problem.  Suppose  you  have  the  m-vector  v  which  is  a  function  of 
the  n-vector  x,  and  further  suppose  that  you  find  it  necessary  to  normalize 
v(x)  to  unit  length,  and  find  the  Jacobian  of  the  result.  That  is,  you  must 
compute  dg/dx  where 


g(x) 


dtf 


v(x) 

VI. 


(55) 


v(x) 

-V(x) 


(56) 


Note  that  g  is  forced  to  be  a  unit  vector,  the  magnitude  of  which  is  inde¬ 
pendent  of  x.  We  desire  the  Jacobian  of  g,  but  it  will  probably  be  virtually 
imposible  to  compute  analytically.  We  are  willing  and  able,  however,  to 
compute  the  Jacobian  of  v.  The  question  is: 


“Is  it  possible  to  obtain  dg/dx  through  a  reasonable  transfor¬ 
mation  of  dv /dx  ?” 


We  can.  indeed,  do  this,  and  this  section  shows  how. 


Define 


and 


Now,  differentiating  Eq.  56.  we  get 

dg , 

dx, 


J  -  - 

dx 

(57 

dx 

( 5^ 

we  get 

(59 

.V2  \  dxj  ‘  dx.,  J 

1  dv ,  1  v,  d.V 

(60 

.V  dx}  .V  .V  dij 

-I  i 


N  \  dxj  dx3  J 


Now  remember;  N  is  the  Euclidean  norm  of  v(x).  That  is. 


so 


N 


U=i 


dN_ 

dij 


U=i 


-1/2 


dvk 

*=1  OX3 


1  A  dvk 

=  v2-v**r 

-¥  *=1  UX1 


A  dvk 

^9kdx} 


k  —  1 


Combining  this  with  Eq.  61  gives 


% 

dx} 


1 

y 


(61) 

(62) 

(63) 

(64) 

(65) 

(66) 


This  is  a  nice,  and  moderately  useful,  result  but  it  admits  of  much  further 
simplification. 

Remember:  from  Eqs.  57  and  58, 


J, 


*j 


dr, 

dxj 


(67) 


anc 


dg, 


'  9:,J  dxj 

We  can  pluck  the  summation  term  out  of  Eq.  66.  and  write  it  as 


(68) 


E  k,  9k 


k~l 


(69) 
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-  EKU 

k= i 

(70) 

n 

C-^ 

c  ^ 

(K  t 

(71) 

Now,  the  term  Jjg  is 
For  the  nonce,  let 

an  n-vector,  of  which  Eq.  71  is 

p  *='  J[g 

the  j-th  component. 

(72) 

Then  Eq.  66  becomes 

dg{  1  (  dvi  { 

dx,~  N\dx,  9,Pl\ 

(73) 

or,  in  terms  of  matrices  instead  of  their  components. 

=  ^{jp-gPr} 

(74) 

=  -g(j,Tg)r} 

(75) 

-  ^{j,.-ggrJ,} 

(76) 

J9  =  ^{l-ggr}jp 

("7) 

Eq.  77  is  the  result  sought.  Is  it  not  truly  fine? 


E  The  Truth  Model 


This  appendix  was  extracted,  virtually  verbatim,  from  the  CLOP  User's  Man¬ 
ual,  and  appears  here  merely  as  reference  material  for  the  reader  who  may 
not  be  familiar  with  how  propagators  of  this  variety  work. 

Since  the  following  two  subsections  were  not  written  expressly  for  this 
report,  this  is  a  little,  minor,  inconsistency  in  the  naming  of  variables.  For 
example,  in  Sec.  E.l,  the  symbol  f  is  used  to  mean  “force  on  the  vehicle”, 
whereas  in  the  rest  of  the  report,  it  means  “nonlinear  differential  equations 
of  motion”.  The  reader  is  beseeched  to  graciously  overlook  these  shortcuts 
that  were  taken  in  the  interest  of  expediency. 


E.l  Differential  Equations  of  Motion 

The  entire  CLOP  program  is  based  on  the  solution  of  the  ubiquitous  vector 
differential  equation 

f  =  ma 


which  can  be  rewritten  as 


where,  of  course. 


In  component  form. 


a  =  r  =  — 
m 


r  = 


dt 2 

m 


/* 

fy 

h  J 


(78) 


(79) 


(80) 


(SI) 


If  the  forces.  /*,  in  this  equation  are  known,  then  it  may  be  integrated  to 
produce  the  position  and  velocity  of  the  vehicle  at  any  time.  At  least,  this  is 
possible  in  principle.  It  is  even  possible,  in  practice,  if  the  forces  are  simple 
enough.  For  instance,  if  the  force  f  results  from  the  gravitational  field  of  a 
spherically  symmetrical  mass  distribution,  then  it  has  the  form 


f  -  -- 


x 

y 


(82) 


1 1 


and  Eq.  81  can  be  integrated  to  produce  the  tidv  elliptical  orbits  of  freshman 
physics  classes.  When  f  becomes  more  complicated  than  this,  as  a  result  of 
modeling  non-symmetric  mass  distributions  and  horrific  atmospheric  drags, 
things  become  intractable  fast.  Then.  Eq.  81  can  no  longer  be  solved  ana¬ 
lytically.  Numerical  methods  must  be  used  to  obtain  approximate  solutions. 
This  is  what  is  done  in  CLOP. 

The  primary  algorithm  used  by  CLOP  to  solve  differential  equations  nu¬ 
merically  is  an  8th  order  Runge-Kutta  process  although  a  4th  order  is  also 
available.  It  is  naturally  suited  to  solving  systems  of  1st  order  simultaneous 
differential  equations;  not  the  2nd  order  set  shown  in  Eq.  81.  To  accommo¬ 
date  R-K,  we  define  a  six  component  System  State  Vector,  u,  as  follows: 


u 


def 


X 

y 

z 

X 


(83) 


The  elements  of  u  are  the  components  of  the  vehicle's  position  and  velocity 
expressed  in  Ec .  Using  this  vector.  Eq.  81  can  be  rewritten  as 


d  u 

Tt 


Au  -  B(ag  -  arf) 


(84) 


where  a?  and  a<*  are  the  accelerations  of  the  vehicle  resulting  from  the  Earth's 
gravitational  field,  and  atmospheric  drag,  respectively.  In  Eq.  84,  A  and  B 
are  constant  matrices.  The  expanded  form  of  Eq.  84,  showing  A  and  B.  is 


i 

■  0 

0 

0 

1 

0 

0  ' 

X 

‘  0 

0 

0  ' 

y 

0 

0 

0 

0 

1 

0 

y 

0 

0 

0 

z 

0 

0 

0 

0 

0 

l 

z 

0 

0 

0 

X 

0 

0 

0 

0 

0 

0 

X 

1 

0 

0 

y 

0 

0 

0 

0 

0 

0 

y 

0 

1 

0 

z 

0 

0 

0 

0 

0 

0 

z 

0 

0 

I 

The  acceleration  forcing  terms  appearing  in  equation  84  are  discussed  in  the 


following  sections. 


I'» 


It  should  be  noted  that  Eq.  84  does  not  contain  ail  the  multitudinous 
terms  which  can  affect  the  acceleration  of  the  vehicle.  As  was  mentioned  in 
this  manual’s  introduction,  the  effects  of  solar  wind.  Earth  albedo,  magnetic 
forces  on  a  charged  satellite,  and  the  gravitational  forces  of  the  Sun,  Moon, 
and  planets  have  not  been  included.  Relativistic  and  astrologic  effects  have 
also  been  ignored. 

E.2  Gravitational  Model 

The  acceleration  a  body  experiences,  as  a  result  of  the  gravitational  field 
in  which  it  is  immersed,  is  a  vector  quantity.  Therefore,  any  mathematical 
model  of  the  gravitational  field  must  be  a  vector  model.  But  the  gravitational 
field  of  the  Earth  (and  virtually  any  other  body  about  which  one  might 
be  interested  in  studying  satellite  orbits)  is  sublimely  complex,  and  simple 
models  of  the  field  will  just  not  do.  The  point  is  this;  the  mathematical 
representation  of  a  complex  vector  field  gets  out  of  hand  quickly. 

Mathematicians,  being  a  basically  lazy  lot,  have  come  up  with  a  simpli¬ 
fication.  It  is  well  known  that  the  gradient  of  a  scalar  function  of  position 
produces  a  vector.9  So,  if  we  can  construct  a  scalar  function,  say  U(r,  A), 
of  position  such  that  its  gradient  matches  the  gravitational  field  of  interest, 
then  we  will  have  only  one  messy  equation  to  tote  about,  and  life  will  be 
much  easier.  Granted,  the  gradient  of  this  function  will  have  to  be  computed 
if  ever  the  function  is  to  be  actually  used  for  anything,  but  that’s  a  mere 
detail  not  pertinent  to  the  theory  of  potential  functions. 

This  scalar  function  of  position  is  commonly  called  a  Potential  Function. 
I  suspect  the  reason  for  this  stems  from  the  fact  that  the  potential  energy 
imparted  to  a  partlcje  by  moving  it  through  a  force  field,  f.  is 

V  -  J  f  «ds  (86) 

and  that  the  gradient  of  this  potential  energy  extracts  the  force;  a  vector 
quantity. 

We  follow  the  same  convention  here,  except  that,  as  the  implementors,  we 
must  compute  the  gradient.  This  is  actually  done  analytically  within  CLOP: 

9This  is  a  matter  of  definition,  actually. 


-10 


numerical  differentiation  is  not  done.  The  gravitational  potential,  V ,  and  its 
gradient  are  treated  in  the  following  sections. 


E.2.1  The  Gravitational  Potential,  f/(r,0,  A) 

The  gravitational  acceleration  term,  ag,  is  derived  as  the  gradient  of  a  scalar 
potential  function,  U .  That  is. 


ag  =  VD' 


In  rectangular  coordinates,  the  gradient  is 


VUC 


dU  ' 
dz 

au 

dy 

au 

di  . 


(87) 


(88) 


Because  of  the  spherical  symmetry  of  the  situation,  however,  the  potential 
function  is  most  easily  expressed  in  terms  of  local  spherical  coordinates,  in 
which  case  the  gradient  becomes 


VUl  = 


au  -] 

dr 


i  au 

rcos[tp)  d\ 


IdU 
r  d<t  J 


(89) 


In  this  form,  r  is  the  radial  distance  from  the  center  of  the  gravitating  body. 
A  is  the  fast  longitude,  and  <j)  is  the  north  geocentric  latitude. 

The  potential  function,  expressed  in  spherical  coordinates  at  thi  point 
whtri  the  vehicle  is  located ,  is 


r  = 


where 


n  =  2 


1  -  Y,  (  J  )  Y1  P^(S(P){Cnm  cos(mA)  -  Snm  sin(mA)} 


m  =0 


(90) 


so  -  sin(c>)  (91) 

The  P™(so)  terms  are  the  Legendre  Associated  Functions  of  degree  n  and 
order  m.  These  functions  are  treated  more  fully  in  Section  E.2.3. 
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Equation  90  is  not  the  result  of  a  late  night  fancy  resulting  from  ov¬ 
erindulgence  at  the  evening’s  supper.  It  is,  rather,  a  solution  to  Laplace's 
equation. 

V2U-—'—  — 

dx2  dy 2  dz2 

in  spherical  coordinates.  This  topic  is  treated  in  excruciating  detail  in  Chap¬ 
ter  II  of  the  book 

"The  Theory  of  Spherical  and  Ellipsoidal  Harmonics” 

E.  W.  Hobson.  Sc.D..  LL.D.,  F.R.S. 

Cambridge,  at  the  University  Press 
1931 

In  essence,  Eq.  92  is  expressed  in  spherical  coordinates,  and  a  solution  to  the 
resultant  partial  differential  equation  is  assumed.  This  assumed  form  leads 
to  several  ordinary  differential  equations  which  must  be  solved.  One  of  the 
differential  equations  is  a  function  of  r  only.  The  other,  known  as  “Legen¬ 
dre's  Equation”,  includes,  as  part  of  its  solution,  the  “Legendre  Polynomials 
and  Associated  Functions”,  P™( sin(d>)).  The  Cnm  and  Snm  terms  are,  es¬ 
sentially,  the  constants  of  integration  of  that  solution.  They  are  determined 
experimentally  by  flying  satellites,  observing  how  they  move,  and  choosing10 
the  coefficients  so  the  predicted  orbits  match  the  observed  ones  as  closely  as 
possible. 

E.2.2  Gradient  of  U(r.  <t>.  A) 

Earlier,  it  was  stated  that  analytical  expressions  for  the  gradient  of  L'(r,0.  A) 
are  coded  into  CLOP.  These  expressions  are  given  here,  purely  for  complete¬ 
ness  of  the  presentation.  In  an  attempt  to  keep  things  under  control,  two 
simplifying  symbols  are  defined.  Let 

rI(m.A)<<=/  {Cnm  cos(  mA )  -  5nm  sin(mA)}  (93) 

and 

f2(m.  A)  de=  {Cnm  sin(rnA)  >nmcos(mA)}  (94) 


°Through  a  highly  sophisticated  Least  Squares  process. 


I 


Then 


dU_ 

dr 

_1 _ au 

r  cos(0)  dX 
1  dU_ 

r  d<j> 


1*2 


n  =  2 

N 

r2  cos (0)  [ 


m =0 


E(-)”E  <(«.,)) 

n=2  '  r  '  m=0 


ficos(<t>) 


’  N  (t  \  n  n  dPm 

E(-) 

Ln=2  V  r  /  m^O  d.sM4>)\ 


(95) 

(96) 


(97) 


E.2.3  Legendre  Polynomials 

The  first  few  Legendre  Polynomials  and  Associated  Functions  will  be  pre¬ 
sented  here  so  the  reader  may  gain  a  feel  for  what  these  creatures  look  like. 
Then,  the  recursion  formulae  CLOP  uses  to  evaluate  these  functions  will  be 
given.  The  symbol  2  will  be  used  to  denote  the  sine  of  the  North  Latitude 
of  the  radius  vector  to  the  vehicle. 


Poi*)  =  1 

P,°(c)  --  2 

*?(*)  =  “(3c:2  -  1) 

P3°(c)  =  ! (5c2  -  3) 

P°(c)  =  ^(3oc4  -  30c2  —  3) 

O 

P}(  =  )  (1  - 

p2'(c)  -  32(1  -  c2)i 


(98) 

(99) 

(100) 

(101) 

(102) 

(103) 
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(104) 


P?(z) 

=  3(1  -  3a) 

(105) 

p*(*) 

=  l(5^2  -  1)(1  -  -2)^ 

(106) 

*?(*) 

=  15z(l  -  z2) 

(107) 

p&) 

=  15(1  -  r2)^ 

(108) 

p}{=) 

=  lz(7z2-m->2)* 

(109) 

p?(z) 

(no) 

p?(=) 

=  105r(l-;2)? 

(111) 

Ptu) 

=  105(1  -  ;2)2 

(112) 

Many  recursion  formulae  exist  for  the  rapid  evaluation  of  these  functions 

and  their  derivatives.  Those  used  by  CLOP  depend  on  the  value  of 
apply  for  those  cases  in  which  n  >  2.  Here  they  are. 

m,  and 

m  =  U 

Pn  = 

( (2n  -  1  )zP°.{  -  (n  -  1  )Pn°.2)/n 

(113) 

dP° 

dz 

dP° 

,  a  rn- 1  nO 

;  dz  -nP-' 

(114) 

m  -  a 

pm  _ 
r  m 

(2m  -  \  )yPZ:\ 

(115) 

dP r 

dz 

/  dPm_.1  -  \ 

,2m  ,T  dz'  '9p-:n 

(1161 

•>U 


Otherwise 


Pn  =  (117) 

d  Pm 

~  =  ((»-m )f?-,  -  nzP?)  is?  (118) 

where 

y  =  \/l  -  ;2  =  cos(d>)  (119) 


F  The  Problem  With  Keplerian  Represen¬ 
tation  of  Orbits 

In  Sec.  2.3.1,  the  state  vector,  u,  was  described  as  consisting  of  the  normal 
Keplerian  orbit  elements;  see  Eq.  1  on  page  12.  Keplerian  orbit  elements 
assume  the  Earth  is  a  perfect  sphere,  and  assume  point-masses  with  their 
associated  uniform  /x/r  potentials.  The  Keplerian  elements  serve  beautifully 
to  describe  the  motion  of  orbiting  objects,  and  the  quantities  are  fairly  easy 
to  compute. 

As  a  first  attempt  at  modeling  a  non-spherical  Earth,  the  general  gravita¬ 
tional  potential,  Eq.  90  page  47,  is  simplified  by  setting  N  =  2,  and  keeping 
only  the  C2o  (  —  J2,  equivalently)  term.  But,  to  solve  this  equation,  numerous 
simplifying  assumptions  must  be  made  to  obtain  any  kind  of  solution  at  all. 
The  result  is  the  “First  Order  Theory”  (FOT)  described  in  Eqs.  2.  3.  and  4. 
The  implication  is  that  these  drift  rates  are  constant.  They  aren't,  especially 
in  the  case  of  du>/dt.  The  theory  also  predicts  that  the  semi-major  axis,  a, 
the  eccentricity,  e.  and  the  orbit's  inclination,  i,  will  be  constant.  This  is 
almost  true.  It  is  close  to  true  if  we  interpret  the  results  as  effects  which  ap¬ 
ply  over  a  single  orbit,  which,  after  all,  are  the  conditions  under  which  they 
were  derived.  The  problem  arises  for  an  AutoNav  system  because  they  do 
not  describe  the  instantaneous  motion  well  at  all.  especially  for  near-circular 
orbits. 

The  FOT  is  useful  for  gaining  broad-brush  insights  into  orbital  motion, 
but  not  for  accurately  predicting  satellite  positions.  This  can  be  seen  in 
Fig.  18.  page  65.  which  shows  how  a  vehicle  flying  according  to  FOT  predic¬ 
tions  departs  from  motion  in  16-th  order  gravitational  field  with  no  air-drag 
or  other  applied  forces.  It  is  seen  that,  for  this  particular  circumstance,  a 
position  error  of  17  km  builds  after  only  one  orbit,  and  things  get  worse  from 
there. 


G  Air  Force  Bases  in  the  U.S 


The  following  is  a  list  of  location  coordinates  for  Air  Force  bases  located  in 
the  political  United  States.  The  original  information  was  obtained  through  a 
literature  search  by  the  BASG  technical  library.  The  list  of  base  coordinates 
is  given  here,  followed  by  the  AWK  program  which  produced  it  from  the 
original  form. 


Air  Fore#  Batts  in  ths  Political  Unitad  Statas 


Lat 

Lon 

Haight 

Location 

(dag) 

(dag) 

(ka) 

35.960 

270.060 

0.077 

!  BLYTHEVILLE 

AK 

34.760 

287.717 

0.094 

!  LITTLE*R0CK 

AK 

32.400 

273.717 

0.067 

!  GOITER 

AL 

32.367 

273.700 

0.061 

!  MAXWELL 

AL 

84.833 

212.900 

0.183 

!  EIELS0N 

AS 

81.260 

210.183 

0.036 

!  ELMEND0RF 

AS 

32.183 

249.117 

0.799 

!  DAVIS-M0NTHAN 

AZ 

33.560 

247.817 

0.332 

!  LUKE 

AZ 

33.360 

248.167 

0.422 

!  WILLIAMS 

AZ 

39.117 

238.633 

0.034 

•  BEALE 

CA 

37.383 

239.433 

0.067 

!  C'  TLE 

CA 

34.900 

242.133 

0.702 

!  El  ■  -ADS 

CA 

34.583 

242.633 

0.876 

•  GEORGE 

CA 

34.087 

241.750 

0.029 

•  LQS*ANGLES 

CA 

33.900 

242.750 

0.466 

!  MARCH 

CA 

38.660 

238.717 

0.029 

!  MATHER 

CA 

34.100 

242.750 

0.352 

!  NORTON 

CA 

38.287 

238.083 

0.019 

!  TRAVIS 

CA 

34.683 

239.517 

0.122 

!  7ANDENBERG 

CA 

38.933 

256.383 

1.910 

!  FALCON 

CO 

39.717 

255.117 

1.846 

!  LOWRT 

CO 

38.817 

255.283 

1.890 

!  PETERSON 

CO 

38.983 

256.133 

2.219 

!  US* AIR*F0RCE* ACADEMY 

CO 

39.117 

284.517 

0.009 

!  DOVER 

DE 

30.483 

273.500 

0.702 

•  EGLIN 

FL 

25.500 

279.600 

0.002 

!  HOMESTEAD 

FL 

30.417 

273.317 

0.011 

!  HULBURT 

FL 

27.850 

277.500 

0.002 

!  MACDILL 

FL 

28.233 

279.400 

0.003 

!  PATRICK 

FL 

30.150 

274.350 

0.005 

!  TYNDALL 

FL 

30.850 

276.750 

0.071 

!  MOODY 

GA 

21.333 

202.100 

0.000 

!  HICX1M 

HA 

21.500 

201.960 

0.258 

!  WHEELER 

HA 

43.050 

244.133 

0.914 

!  MOUITAIN*HOME 

ID 

40.300 

271.850 

0.224 

!  CHJjniTE 

IL 

38.533 

270.133 

0.138 

!  SCOTT 

IL 

37.833 

262.733 

0.418 

!  MCCONNELL 

KS 

32.500 

266.367 

0.051 

!  BARKSDALE 

LA 

31.283 

267.517 

0.027 

!  ENGLAND 

LA 

42.483 

288.717 

0.041 

!  HAISCOM 

MA 

41.850 

289.450 

0.040 

*  QTIS*ANG*BASE 

MA 

42.183 

287 . 433 

0.074 

!  UESTOVER 

MA 

38.800 

283.133 

0.086 

!  A ID REUS 

MD 

48.917 

292.167 

0.230 

!  LORING 

ME 

46.233 

275.533 

0.372 

!  KI*SAVTER 

MI 

42.817 

277.183 

0.178 

!  SELFRIDGE* ANG*BASE 

MI 

44.450 

276.600 

0.193 

!  WURTSMITH 

MI 

38.717 

266.617 

0.265 

!  WHITEMAN 

MO 

33.650 

271.560 

0.065 

.'  COLUMBUS 

MS 

30.417 

271.083 

0.008 

!  XEESLER 

MS 

47.517 

248.800 

1.074 

!  MALMSTROM 

MT 

35.133 

281.017 

0.066 

!  POPE 

NC 

35.367 

282.033 

0.033 

(  SEYMOUR* JOHNSON 

NC 

47.983 

262.967 

0.278 

!  GRAHD* FORKS 

ND 

48.417 

258.667 

0.508 

!  MINOT 

ND 

41.150 

264.050 

0.319 

•  OFFUTT 

NE 

43.067 

289.217 

0.031 

!  PEASE 

NH 

40.017 

285.400 

0.041 

!  MCGUIRE 

NJ 

34.400 

256.733 

1.309 

!  CANNON 

NM 

32.867 

253.900 

1.248 

!  HOLLOMAN 

NM 

35.033 

253.383 

1.631 

!  XIRTLAND 

NM 

36.200 

244.917 

0.570 

!  NELLIS 

NV 

43.217 

284.567 

0.154 

!  GRIFFISS 

NY 

44.867 

286.550 

0.072 

!  PLATTSBURGH 

NY 

39.817 

277.050 

0.227 

!  RICXENBACKER*ANG*BASE 

OH 

39.783 

275.950 

0.251 

!  WRIGHT ‘PATTERSON 

OH 

34.660 

260.683 

0.419 

!  ALTUS 

OK 

35.483 

262.467 

0.393 

!  TINKER 

OK 

36.417 

262.133 

0.398 

!  VANCE 

OK 

32.900 

279.933 

0.014 

!  CHARLESTON 

SC 

33.700 

281.117 

0.008 

!  MYRTLE*BEACH 

SC 

33.987 

279.517 

0.074 

!  SHAW 

SC 

44.167 

256.900 

0.975 

!  ELLSWORTH 

SD 

35.350 

273.800 

0.305 

!  ARNOLD 

TN 

3-1 


30.217 

262.333 

0.165 

!  BERGSTROM 

TI 

29.3S0 

261.567 

0.183 

!  BROOKS 

TI 

32.750 

262.567 

0.198 

!  CARSWELL 

TI 

32.417 

260.200 

0.545 

!  DTESS 

TI 

31.467 

259.550 

0.572 

•  GOODFELLCW 

TI 

29.383 

261.417 

0.210 

!  KELLT 

TI 

29.383 

261.400 

0.227 

•  LACKLAND 

TI 

29.350 

259.167 

0.329 

!  LAUGHLXN 

TI 

33.583 

269.150 

1.017 

!  REESE 

TI 

33.983 

261.483 

0.309 

•  SHEPPARD 

TI 

41.233 

248.050 

1.459 

!  HILL 

0T 

37.083 

283.650 

0.003 

!  LANGLET 

VA 

47.633 

242.367 

0.760 

!  FAIRCHILD 

WA 

47.117 

237.417 

0.098 

•  MCCHORD 

UA 

41.133 

256.183 

1.878 

!  FRANCIS*E» WARREN 

WT 

There  are  89  bases  in  this  fils 

Minimum  latitude:  21.333  degrees 

Maximum  latitude:  84.033  degrees 

Minimum  longitude:  201. 960  degrees 
Maximum  longitude:  292.167  degrees 
Minimum  height :  0 . 000  metres 

Maximum  height:  2218.944  metres 


BEGXB 


BF  ==  0 


END 


{  latmin  =  1.0«10 

latmax  =  -1.0*10 
lonmin  =  1.0*10 

lonmax  =  -1.0*10 
hmin  =  1.0*10 

hmax  -  -1.0*10 


> 


print! ("!\t\t Air  Force 
printK"  !\n”) ; 

Bates 

in  the  Political  United  StateiXn" 

print! ( " ! 

Lat 

Lon 

Height \t  LocationXn") ; 

printK"! 

(deg) 

(deg) 

(km)\n\n") ; 

{  latitude  =  $4  +  18  /  60 . 0 ; 

longitude  =  380.0  -  ($7  +  10  /  80.0); 

height  =  $3  *  3.048*-4; 


latmin  =  latitude  <  latmin  7  latitude 
latmax  =  latitude  >  latmax  ?  latitude 


latmin; 
latmax ; 


lonmin  =  longitude  <  lonmin  7  longitude  :  lonmin; 
lonmax  -  longitude  >  lonmax  7  longitude  :  lonmax; 


hmin  =  height  <  hmin  7  height 

hmax  =  height  >  hmax  7  height 


:  hmin; 
:  hmax ; 


-♦bases; 


printlC7.10.3f  '410.3!  '/.10.3f\t!  ‘/.-24s  7,s\n”, 


} 


latitude , 
longitude, 
height , 

$1,  $2  ); 


{  print! ("\n\n") ; 
print! ("!  There  are  '/.d  bases 
print! (" ! \n") ; 
print! ("!  Minimum  latitude: 
print! ("!  Maximum  latitude: 
printK"!  Minimum  longitude: 
printK"!  Maximum  longitude: 
printl("!  Minimum  height : 
pnntl("!  Maximum  height : 


in  this  !ile\n", 

bases 

); 

'48. 31  degreetXn" , 

latmin 

) 

'48.31  degreetXn", 

latmax 

) 

'48.31  degreetXn", 

lonmin 

) 

'48.31  degreetXn", 

lonmax 

) 

'/.8 . 31  metresXn"  , 

hmin  • 

1000.0) 

‘48.3!  metresXn”, 

hmax  * 

1000.0) 

> 


H  Symbols 

a  Orbit  semi-major  axis 

e  Orbit  eccentricity 

i  Orbit  inclination 

f(u,<)  Nonlinear  differential  equations  of  motion 
F  Linearized  plant  geometry 

g*.  Unit  vector  parallel  to  v*. 

h(u.t)  Nonlinear  measurement  function 
H  Linearized  measurement  geometry 

I  Identity  matrix 

J  Jacobian:  taken  with  respect  to  u 

J2  Second  gravitational  potential  zonal  coefficient 

K  Kalman  gain  matrix 

/x  Earth’s  gravitational  parameter 

v  Orbit  true  anomaly 

lj  Orbit  argument  of  perigee 

fl  Orbit  right  ascension  of  ascending  node 

P  State  covariance  matrix 

Q  Plant  covariance  matrix 

R  Measurement  covariance  matrix 

re  Equatorial  radius  of  the  Earth 

r  Orbit  period 

1 0  Time  at  Epoch 

t  Time 

Ik  Time  at  the  k- th  event 

u  State  vector 

uk  State  vector:  Keplerian  elements 

uc  State  vector:  Cartesian  elements 

Vfc(u,t)  Vector  from  vehicle  to  reference  point  k 
z  Measurement  vector 


NORTH  AMERICA 


EAST  LONGITUDE  (dag) 


Figure  4:  Ground  Track  of  Test  Orbit 
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Figure  7:  Flight  la 


Figure  8:  Flight  lb 
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I 


Figure  9:  Flight  lc 


